library(dplyr)

rm(list = ls())

# load data
load("Brader_et_al2008.RData")

# function for demeaning
demean <- function(x) x - mean(x, na.rm = TRUE)

# function for residualizing intermediate confounders
residualize <- function(formula, df) residuals(lm(formula, df))

# data preprocessing
Brader2 <- Brader %>%
  select(immigr, emo, p_harm, tone_eth, ppage, ppeducat, ppgender, ppincimp) %>% na.omit() %>%
  mutate(immigr = 4 - immigr,
         hs = (ppeducat == "high school"),
         sc = (ppeducat == "some college"),
         ba = (ppeducat == "bachelor's degree or higher"),
         female = (ppgender == "female")) %>%
  mutate_at(vars(emo, p_harm, ppage, female, hs, sc, ba, ppincimp), demean) %>%
  mutate(., p_harm_res = residualize(p_harm ~ ppage + female + hs + sc + ba + ppincimp + tone_eth, .))

# total effect model
total_mod <-  lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth,
                 data = Brader2)

# naive estimate of direct effect

naive_mod <- lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth + 
                  emo * (ppage + female + hs + sc + ba + ppincimp + tone_eth), data = Brader2)

summary(naive_mod)

# sequential g

ptbias_mod <- lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm + 
                   emo * (ppage + female + hs + sc + ba + ppincimp + tone_eth),
                 data = Brader2)

m_indices <- c("emo", "ppage:emo", "female:emo", "hs:emo", "sc:emo", "ba:emo", "ppincimp:emo", "tone_eth:emo")

y_demed <- model.response(model.frame(ptbias_mod)) -
  model.matrix(ptbias_mod)[, m_indices] %*% coef(ptbias_mod)[m_indices]

direct_mod <- lm(y_demed ~ ppage + female + hs + sc + ba + ppincimp + tone_eth, data = Brader2)

# rwr without intermediate interactions
rwr1_mod <-  lm(immigr ~  ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res +
                  emo * ( ppage + female + hs + sc + ba + ppincimp + tone_eth),
                data = Brader2)

# rwr with intermediate interactions
rwr2_mod <-  lm(immigr ~  ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res +
                  emo * (ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res),
                data = Brader2)

# bootstrap 

nboots <- 500

total_hold <- matrix(NA, nrow = length(coef(total_mod)), ncol = nboots)
naive_hold <- matrix(NA, nrow = length(coef(naive_mod)), ncol = nboots)
direct_hold <- matrix(NA, nrow = length(coef(direct_mod)), ncol = nboots)
rwr1_hold <- matrix(NA, nrow = length(coef(rwr1_mod)), ncol = nboots)
rwr2_hold <- matrix(NA, nrow = length(coef(rwr2_mod)), ncol = nboots)

set.seed(02138)

for (i in 1:nboots) {
  
  star <- sample(1:nrow(Brader2), replace = TRUE)
  Brader2_star <- Brader2[star, ]
  
  Brader2_star <- Brader2_star %>% tbl_df() %>%
    mutate(., p_harm_res = residualize(p_harm ~ ppage + female + hs + sc + ba + ppincimp + tone_eth, .))
  
  total_star <-  lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth, data = Brader2_star)
  
  naive_star <- lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth + 
                     emo * ( ppage + female + hs + sc + ba + ppincimp + tone_eth), data = Brader2_star)
  
  ptbias_star <- lm(immigr ~ ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm + 
                      emo * ( ppage + female + hs + sc + ba + ppincimp + tone_eth),
                    data = Brader2_star)
  
  m_indices <- c("emo", "ppage:emo", "female:emo", "hs:emo", "sc:emo", "ba:emo", "ppincimp:emo", "tone_eth:emo")
  
  y_demed <- model.response(model.frame(ptbias_star)) -
    model.matrix(ptbias_star)[, m_indices] %*% coef(ptbias_star)[m_indices]
  
  direct_star <- lm(y_demed ~ ppage + female + hs + sc + ba + ppincimp + tone_eth, data = Brader2_star)
  
  rwr1_star <-  lm(immigr ~  ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res +
                     emo * ( ppage + female + hs + sc + ba + ppincimp + tone_eth),
                   data = Brader2_star)
  
  rwr2_star <-  lm(immigr ~  ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res +
                     emo * (ppage + female + hs + sc + ba + ppincimp + tone_eth + p_harm_res),
                   data = Brader2_star)
  
  total_hold[, i] <- coef(total_star)
  naive_hold[, i] <- coef(naive_star)
  direct_hold[, i] <- coef(direct_star)
  rwr1_hold[, i] <- coef(rwr1_star)
  rwr2_hold[, i] <- coef(rwr2_star)
}

rownames(total_hold) <- names(coef(total_mod))
rownames(naive_hold) <- names(coef(naive_mod))
rownames(direct_hold) <- names(coef(direct_mod))
rownames(rwr1_hold) <- names(coef(rwr1_mod))
rownames(rwr2_hold) <- names(coef(rwr2_mod))

out_coefs <- c("(Intercept)", "tone_eth", "emo", "tone_eth:emo", "p_harm_res", "p_harm_res:emo")

total_est <- coef(total_mod)[out_coefs]
naive_est <- coef(naive_mod)[out_coefs]
direct_est <- coef(direct_mod)[out_coefs]
rwr1_est <- coef(rwr1_mod)[out_coefs]
rwr2_est <- coef(rwr2_mod)[out_coefs]

total_se <- apply(total_hold, 1, sd)[out_coefs]
naive_se <- apply(naive_hold, 1, sd)[out_coefs]
direct_se <- apply(direct_hold, 1, sd)[out_coefs]
rwr1_se <- apply(rwr1_hold, 1, sd)[out_coefs]
rwr2_se <- apply(rwr2_hold, 1, sd)[out_coefs]

table1_est <- cbind(total_est, naive_est, direct_est, rwr1_est, rwr2_est) %>%
  `row.names<-`(out_coefs)
table1_se <- cbind(total_se, naive_se, direct_se, rwr1_se, rwr2_se) %>%
  `row.names<-`(out_coefs)

est_indices <- 2 * (1:nrow(table1_est)) - 1
se_indices <- 2 * (1:nrow(table1_se))

table1 <- matrix(NA, nrow = 2 * nrow(table1_est), ncol(table1_est))
table1[est_indices, ] <- table1_est[, ]
table1[se_indices, ] <- table1_se[, ]

rownames(table1) <- rep(NA, 12)
rownames(table1)[est_indices] <- out_coefs
rownames(table1)[se_indices] <- "S.E."
write.csv(table1, file = "table1.csv")
